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Abstract. The problem of characterizing complexity of quantum dynamics - in 
particular of locally interacting chains of quantum particles - will be reviewed and 
discussed from several different perspectives: (i) stability of motion against external 
perturbations and decoherence, (ii) efficiency of quantum simulation in terms of 
classical computation and entanglement production in operator spaces, (iii) quantum 
transport, relaxation to equilibrium and quantum mixing, and (iv) computation of 
quantum dynamical entropies. Discussions of all these criteria will be confronted 
c"j . with the established criteria of integrability or quantum chaos, and sometimes quite 

surprising conclusions are found. Some conjectures and interesting open problems in 
q-( ergodic theory of the quantum many problem are suggested. 
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1. Introduction 



This article will be about discussing the possibilities to characterize randomness and 
dynamical complexity in quantum mechanics and relating this issue to the questions of 
non-equilibrium statistical mechanics. We shall try to illustrate, mainly by presenting 
various numerical examples, a possible 'cyclist approach' [1] towards the quantum many- 
body problem which is inspired by experiences gained in studies of quantum and classical 
chaos of one or few particles (see e.g. [21 [1]). 

Solving the many body problem in quantum mechanics presents a major challenge 
from its very beginnings. And along this way, many ingenious important analytical and 
efficient numerical methods of solution have been developed, for example Bethe ansatz 
[3], later generalized and interpreted as the quantum inverse scattering problem [I], 
real space renormalization group methods [5], quantum Monte Carlo [6], density matrix 
renormalization group (DMRG) [7], and more recently various quantum information- 
theoretic based time-dependent DMRG (tDMRG) [8j [9j [10]. The ultimate aim of any 
of these methods is to efficiently find analytical or numerical approximation to the 
solution for some of the physical observables in the quantum many body problem, 
however many methods work only under some specific conditions which are not always 
well understood. For example, Bethe ansatz and quantum inverse scattering work only 
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for a small subset of problems which are completely integrable, and which may often 
have very non-generic non-equilibrium properties, such as e.g. ballistic transport at 
high temperatures [11]. Quantum Monte Carlo techniques represent a very successful 
set of numerical methods which can yield thermodynamic equilibrium averages for 
generic (non-integrable) systems, however they are practically useless for computing non- 
equilibrium, or time dependent quantities, such as transport coefficients. Traditional 
DMRG method |7J is provably very successful for computing accurate ground state 
expectation values of almost any physical observables of one- dimensional interacting 
systems. And tDMRG [3J [10] promised to extend the success of DMRG to time 
dependent physics. Indeed, the first numerical experiments looked very promising, but 
after a closer look one may realize that they have all been applied to a rather special 
subset of interacting systems and to a rather special subset of initial states. For generic 
(non-integrable) interacting systems or for sufficiently complex initial states, tDMRG 
should fail to provide an efficient computation as shall be discussed below. In our view 
this is to be expected, and represents an intrinsic characteristic of quantum complexity 
of such system and should correspond to many body extensions of the phenomena 
of quantum chaos p] where spectra and eigenvector coefficients can be described by 
statistical ensembles of random matrices. In the past two, almost three decades there 
has been a lot of activity in the so called field of quantum chaos, or quantum chaology 
[12] , where people were trying to understand the essential and significant features of 
quantum systems which behave chaotically in the classical limit. Classical chaos can be 
defined in terms of positive (algorithmic, Kolmogorov-Chaitin) complexity of systems' 
orbits. Still, the question whether such definition of complexity can be in an intuitively 
sensible way translated to quantum mechanics failed to be answered in spite of many 
efforts. It seems that quantum systems of finite (one or few body) chaotic systems 
are generically more robust against imperfections [13] than their classical counterparts, 
which is consistent with a simple illustration based on wave-stability of unitary quantum 
time evolution [14|. 

Nevertheless, it has been suggested that exponential sensitivity to initial conditions, 
the essential characteristics of classically chaotic systems, has many fingerprints in 
quantum mechanics. For example, in one of the pioneering works on quantum Loschmidt 
echoes (or fidelity decay), Jalabert and Pastawski [15] have shown that in the semi- 
classical regime, decay of system's sensitivity to external perturbation, as defined 
by fidelity, is exponential with the rate which can be related to classical Lyapunov 
exponents. However, this is only true in sufficiently semi-classical regime, where effective 
Planck constant is smaller than effective strength of perturbation, and where fidelity can 
be essentially explained classically |16j . On the other hand, in purely quantum regimes, 
quantum fidelity decays in completely different manner than the classical fidelity, and 
quite surprisingly displays slower decay for systems with stronger decay of temporal 
correlations [13]. 

Throughout this paper we shall discuss various ergodic properties of a simple toy 
model of non-integrable interacting quantum many body system, namely kicked Ising 
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(KI) chain [17], and its time- independent version, which (both) undergo a transition 
from integrable to quantum chaotic regimes when the direction of the (kicking) magnetic 
field is changed. By simulating the dynamics of the model in terms of tDMRG we find 
that it can be performed efficiently only in the integrable regimes. Being interested in 
statistical mechanics of the model we shall describe numerical experiments addressing 
the questions on the relationship between the onset of quantum chaos in KI model and 
its quantum mixing property, namely the nature of relaxation to equilibrium, and the 
properties of non-equilibrium steady state. On one hand we argue that the regime of 
quantum chaos essentially corresponds to the regime of quantum ergodicity and quantum 
mixing where diffusive transport laws are valid. On the other hand, we conjecture that 
the transition to non-ergodic regime may occur before the system parameters reach 
the integrable point (even in thermodynamic limit), and that non-ergodic to ergodic 
transition can be characterized with order parameters which change discontinuously 
at a critical value of system parameter. We argue that the process of relaxation to 
equilibrium in quantum chaotic (mixing) case can be described in terms of quantum 
analogues of Ruelle resonances [18]. Furthermore, we conjecture that the eigenvectors 
of this process possess a certain scale invariance which can be described by simple 
power laws. We also discuss a possibility of numerical calculation of quantum dynamical 
entropies [T9J [20j EI] in a non-trivial setting of KI model, and find, quite remarkably, 
that positivity of such dynamical entropies does not correspond to any other measures 
of quantum chaos, namely quantum dynamical entropies appear to be positive even in 
the integrable regimes. 

About two thirds of the material presented in this paper constitutes a review of 
a selection of recent results related to quantum chaos in many-body systems, with a 
flavor of quantum information. However, about one third of the material, constituting 
a major part of section is new and original and has not yet been published before. 
The paper is organized as follows. In short section [2] we review basic definitions of 
complete integrability and quantum chaos, in particular in the context of many-body 
systems. In section |3] we introduce KI model which will serve us as a very convenient 
and efficient toy model to numerically demonstrate all the phenomena discussed in 
this paper. In section H] we discuss the robustness of quantum systems to external 
perturbations and decoherence in open quantum systems, mainly as characterized by 
quantum Loschmidt echoes and entanglement between the system and the environment. 
In section [5] we discuss the time efficiency of best (known) classical simulation, namely 
using tDMRG, of locally interacting ID quantum systems and its possible dependence 
on integrability of the model. In section [6] we relate standard criteria of quantum chaos 
to normal (diffusive) versus anomalous (ballistic) transport and discuss a simulation 
of a quantum heat current in non-equilibrium steady state. In section [7] we discuss a 
problem of quantum relaxation to equilibrium, i.e. the quantum mixing property, and 
quantitative characterization of quantum dynamical complexity. This section contains 
a large portion of original intriguing numerical results which support few interesting 
and perhaps surprising conjectures. In section [8] we conclude and discuss some open 
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problems. 

2. Integrability versus chaos 

The central issue of this paper is to verify and demonstrate to what extend the complete 
integrability affects non-equilibrium properties of quantum many-body systems and their 
dynamical complexities, and conversely, whether (strong) non-integrability generically 
coincides with the established criteria of quantum chaos. 

Let us start by giving some established definitions (see e.g. [H El E21 H] ) of the basic 
terms needed to understand the issue. 

A classical Hamiltonian many-body system with L degrees of freedom is completely 
integrable, if there exist L functionally independent global smooth phase space functions 
(integrals of motion) which are mutually in involution, i.e. all Poisson brackets among 
them vanish. In such a case global canonical transformation to canonical action- angle 
variables can be constructed, and dynamics can be explicitly solved - at least in principle. 
For locally interacting infinite systems (L — > oo) analytic methods for an explicit 
construction of integrals of motion and canonical action-angle variables are known which 
usually go under a common name of inverse scattering method. Explicit solution by a 
mapping to an iso-spectral Schrodinger problem in terms of inverse scattering technique 
is usually understood as a definition of complete integrability in such context. 

Definition of quantum complete integrability is less unique. Nevertheless, algebraic, 
non-commutative versions of the inverse scattering technique exist and can be applied 
to some quantum interacting lattices in one-dimension, generalizing the famous Bethe 
ansatz solution of the Heisenberg spin 1/2 chain, and this is used as the most general 
known definition of quantum complete integrability. Other versions of integrability 
for finite L quantum systems have been proposed but they do not correspond to the 
integrability of the underlying classical limit, if the latter exists, so these ideas will not 
be considered in this paper. 

It has not been generally accepted yet, though demonstrated in many occasions, 
that completely integrable quantum systems constitute only a small subset of 
physical models and posses many exceptional (non-generic) non-equilibrium dynamical 
properties, like for example anomalous transport at finite temperatures (see e.g. [11].) 

On the other extreme of ergodic hierarchy we have chaotic systems. In classical 
Hamiltonian dynamics of few particles, chaos is best defined in terms of positive 
algorithmic (Kolmogorov) complexity of systems' trajectories, or equivalently, by 
exponential sensitivity to initial conditions. However, bounded quantum systems of 
finite number of particles cannot be dynamically complex as their excitation spectrum 
is discrete, and hence the evolution is necessarily quasi-periodic (or almost periodic). 
Still, quite surprisingly, even for such systems certain dynamical properties are random 
and universal, if the underlying classical limit is sufficiently strongly chaotic [12], [2] . But 
genuine dynamical complexity may emerge in thermodynamic limit. However, there 
is still no completely satisfactory definition of dynamical chaos for infinite quantum 
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systems. The commonly used working definition is the reference with the random matrix 
theory [2], namely the many-body quantum system is said to be quantum chaotic if its 
excitation spectrum or some other dynamical properties can be (on certain energy, or 
time scale) well described by ensembles of random Hermitian matrices with appropriate 
symmetry properties. 

3. Toy model 

Throughout this paper we shall use, either for illustration of theoretical results, or 
for numerical experimental studies, the following ID locally interacting quantum lattice 
system, namely a chain of L qubits, or spin 1/2 particles, coupled with nearest neighbour 
Ising interaction and periodically kicked with spatially homogeneous, but arbitrarily 
oriented magnetic field. In a suitable dimensionless units our model can be written in 
terms of a three-parameter periodic time-dependent Hamiltonian 

L-l 

H(J, h x , h z - t) = Y, {JvJvJ+i + S P (t)(h x a] + Mf)} • (1) 

3=0 

d~ p = J2mez — ^) is a unit-periodic Dirac delta and crj ,7,z are the usual Pauli spin 
variables on a finite lattice j ; € = {0, 1, ... L — 1}, satisfying the commutation (Lie) 
algebra [c™, erf] = ^ 2i£ a p 1 8jkcrj ■ Sometimes it will turn fruitful if we extend the set of 
Pauli matrices by identity matrix and assign them a numerical superscript a", a G Z4, 
namely cr° = 1, crj = aj, a? = crj, cr| = crj. The finite chain will often be treated with 
periodic boundary conditions, o~2 = <Tq , and sometimes the thermodynamic limit (TL) 
L — ► 00 will be considered, in particular when we shall study the statistical mechanics 
of (JTj) in section [71 

Although kicked Hamiltonian one-particle models have been very popular in the 
field of nonlinear dynamics and quantum chaos for decades, for example the Chirikov's 
kicked rotator model [23] , the use of kicked systems in quantum many body physics has 
been so far very limited. If one is not only interested in zero temperature (ground state) 
or low temperature physics, then as we shall try to demonstrate in this review, kicked 
many-body models like ([1]) provide simpler and clearer phenomenological picture of 
global dynamics than time-independent models. The main reason is that since energy 
is not a conserved quantity, the entire Hilbert space of many-body configurations is 
accessible for non-integrable dynamics, and the notions of quantum ergodicity and 
mixing (see e.g. Ref.[2l] for definitions and further references) can be defined more 
clearly than in the time- independent, autonomous setting. Perhaps the first kicked 
interacting quantum lattice has been introduced in Ref.[25]. Even if in traditional solid 
state physics such kicked dynamics would represent un-physically strong excitations, 
one has to realize that kicked quantum chains could be attractive options as benchmark 
models of quantum state manipulation and quantum computation in optical lattices. 

The Kicked Ising (KI) model ([1]) has been defined for the first time in Ref.|17j. 
generalizing the integrable KI model with transverse field introduced in Ref . [26] . Clearly, 
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as shown there [26J for the case of transverse field, h x = 0, the time-dependent model 
(TjQ) can be considered completely integrable since it can be solved explicitly, for example 
by Wigner-Jordan-Bogoliubov transformation to non-interacting spinless fermions, and 
a large class of its time correlation functions can be calculated explicitly. In addition, an 
infinite sequence of local conserved quantities (integrals of motion) can be constructed 
in such a case. There is another, trivial completely integrable regime of KI model, 
namely the case of longitudinal field, h z — 0. Yet another, more non-trivial completely 
integrable regime of KI model is found when the magnitude of dimensionless field is a 
multiple of 7r/2, namely h = \/h x + h\ = nir/2,n G Z, since then the magnetic kick 
can be considered as a multiple of n rotation, and generated by a slightly obscure set of 
non-interacting spinless fermions. 

However, in a general case of titled magnetic field when both components h x 
and h z are non- vanishing, and 2h/n is non-integer, the model is non-integrable, and 
is conjectured to be not amenable to exact analytical methods. As discussed in the 
following sections, non-integrable KI model can display a variety of regimes according 
to the criteria of quantum chaos, quantum ergodicity and quantum mixing. In fact, 
recently the spectral statistics of KI model in strongly non-integrable regime has been 
studied in detail and random matrix behaviour has been clearly confirmed at short 
energy ranges [27]. 

Due to kicked nature of interaction, the evolution propagator of KI model for one 
unit period of time (the so-called Floquet operator), starting just before the kick, can 
be constructed explicitly in terms of a time-ordered product 

U(J, h x , h z ) = T exp f-i J dtH(J, h x , h z - 1 - 0)^ (2) 

= exp(-iJ^aJa^ fl )exp(-i^(/i x( rJ + /i z (T J z )) (3) 

j j 

The last line suggests a simple efficient quantum protocol to simulate KI model in 
terms of simple 1-qubit U'(h x , h z ) = exp(— i(/i x cr x + h z a z )) and 2-qubit U"(J) = 
exp(— iJo" x ® cr x ) quantum gates. If we write a compact Kicked Ising 2-qubit gate 

W(J, h K , h z ) = U"(J) ■ (U'(h x , h z ) ® 1) (5) 

and introduce a left-to-right ordered operator product, namely Y[j Aj = • • • A 1 A 2 A 3 ■ ■ •, 
then we can write KI protocol as a simple string of W— gates 

+ 

u(j, h x , h z ) = n w jJ+1 (j, h x , h z ). (6) 

3 

The KI model has a better known autonomous limit, namely time-independent 
Ising chain in a tilted magnetic field 

H'(J, h x , h z ) = lim -H(rJ,Th x ,Th z ;t) = ^{Jrf(J* +l + h x a] + h z aj), (7) 

j 
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which is again non-integrable unless the field is transverse h x = 0, or longitudinal h z = 0. 

4. Decoherence and fidelity 

4-1. Loschmidt echoes 

One of the central questions about the dynamics of complex quantum systems is their 
robustness against small imperfections in the Hamiltonian. While it is clear that due 
to unitarity the quantum evolution is stable against variation of initial states [13], it is 
not so clear how stable it is against variation of the Hamiltonian, either being static or 
time-dependent - perhaps even noisy. 

Let us write the Hamiltonian as H$ = Hq + 5V, where Ho is the unperturbed 
Hamiltonian, V is a Hermitian perturbation operator and 5 is a small perturbation 
parameter. Peres [28] proposed the following measure of stability of quantum evolution: 
Let us start from some fixed initial state and write the time evolutions of this state 
generated with the unperturbed and perturbed Hamiltonian, as \ipo(t)) = Uo(t)\ip) and 
\if>5(t)) = Us(t)\ip), respectively. Then the stability is characterized by fidelity, i.e. the 
squared overlap between these two states 



There are two alternative interpretations of quantum fidelity: (i) One can interpret 
(JE]) as a quantum Loschmidt echo, namely the probability that the initial state 
which is propagated forward with perturbed evolution Us(t), and after that propagated 
backwards in time with unperturbed evolution U^t) = U (—t), is again measured in 
the same (initial) state. Alternatively (ii) is just the square modulus of expectation 
value of the unitary echo operator U (—t)Us(t) which is the quantum propagator in the 
interaction picture. 

There have been three main theoretical approaches to understanding of fidelity 
decay in quantum dynamical systems: 

4-1.1. Semi- classical approach. In a seminal work Jalabert and Pastawski derived 
quantum Loschmidt echo for systems which posses well defined classical limit. They 
have shown that under certain conditions, namely that the perturbation strength is 
large enough - typically larger than appropriately scaled Planck constant, and that 
initial states have certain classical interpretations - like coherent states, position states, 
etc, the quantum Loschmidt echo decays exponentially 



with the rate A which is perturbation independent and typically equals the Lyapunov 
exponent of the underlying classical dynamics. More recently a completely classical 
interpretation of this so-called Lyapunov decay has been given [16] in terms of the 
classical Loschmidt echo and a theory for exponents A has been developed in terms of 
the full spectrum of Lyapunov exponents for classical dynamical systems with few [16] 
or many [29J degrees of freedom. 



F(t) = \(Mt)\Mt))\ 2 = mulimmw- 



(8) 




(9) 
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4-1.2. Time dependent perturbation theory and linear response approximation. In the 
opposite regime, where the scaled Planck constant is bigger than the perturbation 
strength 5, one can use time-dependent quantum perturbation theory to second order 
to derive a simple linear response formula for fidelity decay [TTJ, [13] 

F(t) = 1 - ^ J* dt' J dt"C(t', t") + 0{{5/hf) (10) 

in terms of 2-point time correlation function of the perturbation 

C(t\t") = (V(t')V(t")} - (V(t'))(V(t")} (11) 

where V(t) = U (—t)VU (t) is the perturbation operator in the interaction picture and 
(.) = is an expectation value in the initial state. 

From this formula - which can be viewed as a kind of Kubo-like linear response 
theory for dissipation of quantum information, an interesting conclusion can be drawn: 
Fidelity decays faster for systems with slower decay of temporal correlations, or 
alternatively phrased, quantum system is more robust against external perturbations 
if it relaxes to equilibrium faster. 

One can actually go beyond the second order perturbation theory, and sum up 
the Born series for fidelity to all orders in many specific situations |13J . Since we are 
here mainly interested in qubit (spin 1/2) chains, we shall only review specific results 
for Kicked Ising chain [T7j. We shall discuss three different specific values of system 
parameters, in all three we take J = 1.0, h z = 1.4: (a) Integrable regime of transverse 
field h x = 0, (b) weakly non-integrable regime with h x = 0.4, and (c) strongly non- 
integrable regime with h x = 1.4. All the results, correlation functions and fidelity, 
are for random initial states, which can be interpreted as pure states of maximum 
quantum information, and averaging over random states is equivalent to a tracial state 
(.) = 2~ L tr(.). We consider the evolution operator ([3]) and perturb it by changing the 
magnetic field such that the perturbation is generated by the transverse component of 
the magnetization M = J^. cr|, namely Ug(t) = [U(J, h x , h z ) exp(— iSM)] 1 . 

In the integrable case, or in general, in non-ergodic case, where the correlation 
function C(t',t") = (M(t')M(t")) (note that (M(t)) = 0) approaches a non-vanishing 
plateau value as \t' — t"\ — > oo, namely time averaged fluctuation of magnetization 
D M = lim t _ +00 (l/t 2 ) J * dt' f dt"C(t', t") is non-vanishing D M ^ 0, one can sum up Born 
series to all orders, yielding a Gaussian decay of fidelity 

F(t) = exp(-(t/r ne ) 2 ), r nc = D" 1 / 2 *" 1 . (12) 

The only assumptions are that t is long enough for the time average in the definition of 
Dm to converge, and short enough that fidelity is still above the finite size plateau value 
F* ~ 1/2 L . Note that Dm can be considered as an analog of Drude weight or charge 
stiffness in the linear response solid state transport theory (see e.g. [TT]). In figs.( |T|2l 
(a,b) one can observe the correlation plateaus of correlation functions for integrable and 
weakly non-integrable cases and respective good agreement with a Gaussian decay of 
fidelity (112p . Note that in the integrable case the correlation function and the plateau 
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Figure 1. Correlation decay for three cases of finite KI model: (a) integrable h x = 0, 
(b) intermediate h x = 0.4, and (c) ergodic h x = 1.4, whereas J = 1.0, h z = 1.4, and for 
different sizes L = 20, 16, 12 [solid-dotted connected curves, almost indistinguishable 
in (a,b)]. Circles in (a) show exact result for L — oo. 



Dm have been calculated also analytically [26J. Note that for comparing different lattice 
sizes L a size-scaled value of the perturbation strength 5' = 5 \J L/L , with L = 24, has 
been fixed rather than 6 itself. 

One the other hand, for sufficiently strong integrability breaking, say in case (c) of 
KI model, the correlation function C{t' — t") = C(t', t") decays to zero in TL, which can 
be interpreted as quantum mixing behviour. It has been found that quantum mixing 
behaviour typically corresponds to random matrix (quasi)energy level statistics, see e.g. 
Rcf. [24J and the subsequent sections of the present paper. We shall call such behaviour 
the regime of quantum chaos. Here again, Born series for fidelity can be summed up to 
all orders, yielding an exponential decay [TT| 

F(t) =exp(-t/r m ), r m = l/(2a5 2 ) (13) 
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Figure 2. Average fidelity amplitude (absolute value of it) for three cases of finite 
KI: (a) integrable h K = 0, (b) intermediate h x = 0.4, and (c) ergodic h x — 1.4, 
whereas J = 1.0, h z = 1.4, and for different sizes L = 20, 16, 12, and different scaled 
perturbations 6'. Chain curves give theoretical predictions [T71 113j . 



where a = J C(t')dt' is a transport coefficient. The assumptions for validity of ([TBI 
are that t is larger than a characteristic mixing time scale on which C(t) decays and 
short enough so that finite size effects in quantum correlation function C(t) does not yet 
affect the transport coefficient a. This regime of fidelity decay is sometimes referred to 
as the Fermi Golden Rule regime. Again, as demonstrated in figs. ( fi"P|) c the agreement 
of numerical data for KI model with the theory is excellent. 

Note that decay time of fidelity scales as oc 5~ 2 for ergodic and mixing dynamics, 
and as oc <5 _1 for non-ergodic dynamics, making the latter much more sensitive to small 
perturbations. 

4-1.3. Random matrix theory (RMT). Complex quantum systems can often be well 
described by statistical ensembles of Hamiltonians (30] . Assuming that H and V both 
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belong to canonical ensembles of Gaussian random matrices, one can evaluate ensemble 
averages of the fidelity and relate C(t) to the spectral form factor of the underlying 
random matrix ensemble. The perturbative (linear response) theory has been developed 
in Ref. [31], see also Ref. j32j[33]. Furthermore a non-perturbative (super-symmetric) 
averaging has been successfully applied to obtain exact expressions for average fidelity 
amplitude in the most interesting cases [31] . RMT theory of fidelity has been successfully 
applied to chaotic quantum systems and even to several experimental situations |35j . 

Following more applied philosophy, groups from Toulouse and Como have performed 
a series of numerical experiments [36] analyzing the robustness of several reasonable 
models of quantum computer hardware under small imperfections, being either due to 
(static) unwanted inter-qubit coupling or due to stochastic (noisy) unwanted coupling to 
the environment, when performing quantum algorithms simulating various toy models 
of classical and quantum single-particle chaos [33], like for example quantum kicked 
rotator [37J. The numerical results are in line with a general linear response theory, 
stating that static perturbations are in general more dangerous than noisy ones. 

4-2. Decoherence and entanglement between weakly coupled systems 

Similar thinking as in the previous subsection can be applied to perhaps even more 
fundamental problem of quantum physics, to the problem of decoherence [3H]. Here 
we shall limit ourselves to an abstract unitary model of decoherence, where we treat a 
complete unitary evolution over two subsystems, a central system C, and an environment 
E, and then address a relevant information about the central system (i.e. the part of 
the system which is of physical interest) by partially tracing over the environment. 

Such a discussion can indeed be followed with a close link to the problem of 
Loschmidt echoes, by writing the total Hamiltonian H$ = H + 5V as an ideal 
(unperturbed) separable evolution H = He <S> 1e + lc ® He perturbed by a small 
coupling V between the system and the environment. We shall also assume that we 
start from initial pure state which is a product state = \ipc) ® IV'e)- We are interested 
in the properties of a generally mixed state of the central system obtained by partial 
tracing over the environment 

p c (t) =tr E [U 5 (-tM)(i;\Us(t)} (14) 

Then, under an ideal evolution, the state of the system remains a product state at all 
times, so the state of the central system remains pure and there is no decoherence. 
Decoherence is usually characterized in terms of decaying off-diagonal matrix elements 
of pc in a suitable pointer state basis, for example in the eigenbasis of the central system 
Hamiltonian He- In fact for certain special forms of the perturbation operator V, the 
magnitudes of off-diagonal matrix elements of pc can be literally written as fidelities, or 
quantum Loschmidt echoes, in the evolution of the environment perturbed by system- 
environment coupling |39j. In such framework, the Lyapunov regime of fidelity decay 
corresponds exactly to the Lyapunov growth of decoherence discussed by Zurek and Paz 
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However, another interesting indicator of decoherence is the growth of bi-partite 
entanglement between the system and the environment. Perhaps this notion is even 
more general since it does not depend on a particular choice of the pointer state basis. 
Since the state of the universe (central system + environment) is always pure, the 
characterization of entanglement is easy, say in terms of Von Neuman entropy of pc, 
S(t) = — tr c p c logpc) or linear entropy 5 2 (t) = — logtr c p^(t) which is a negative 
logarithm of purity P(t) = trp^(t). Usually, the quantities S(t) and S 2 (t) essentially 
give equivalent results - namely they can both be interpreted as the logarithm of an 
effective rank of the state pc, but the linear entropy, or purity, is more amenable to 
analytical calculations. 

Again considering chaotic models for the central system and the environment, Miller 
and Sarkar [41] (see also [42J) have been able to observe the 'Lyapunov regime' of 
entanglement growth, namely S(t) is for sufficiently strong perturbations 5 found to 
grow linearly S(t) = ht with the rate h which is perturbation independent and given 
by (the sum of positive) classical Lyapunov exponents of the unperturbed, separated 
(sub) systems. 

On the other hand, in the purely quantum regime of small perturbation 5, 
again typically smaller than effective Planck constant, one can use time-dependent 
perturbation theory and derive perturbation dependent entanglement entropy [43"t 144]. 
namely the purity can be explicitly expressed as P(i) ~ 1 — 5 2 C(t) where C(t) is a 
particular integrated correlation function of the perturbation. In this regime we again 
find that quantum systems, and quantum environments, which display faster relaxation 
to equilibrium, are more robust against decoherence due to typical couplings. For studies 
of bi-partite entanglement, in particular in KI model system, see Refs. [45, 146] . 

There exists even closer connection to fidelity theory, namely one can prove a general 
inequality [431 IS] stating that purity is always bounded by the square of fidelity 

(F(t)Y < P(t). (15) 

In other words: log(l/F(t) 2 ) gives an upper bound for the growth of the linear 
entanglement entropy, 5" 2 (t) < — 21ogF(t). 

In a slightly different context, bit-wise entanglement between a pair of qubits of a 
quantum register representing a time dependent quantum state, where the rest of the 
register is considered as an environment, has been demonstrated to be an indicator of 
quantum chaos [48] and even a signature of classical chaos [49] . 

5. Efficiency of classical simulations of quantum systems 

In the theory of classical dynamical systems there is a fundamental difference between 
integrable and chaotic systems as outlined in section [3 Chaotic systems, having positive 
algorithmic complexity, unlike the integrable ones, cannot be simulated for arbitrary 
times with a finite amount of information about their initial states. Computational 
complexity of individual chaotic trajectories is linear in time, however, if one wants 
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to describe statistical states (phase space distributions) or observables of chaotic 
classical systems, up to time t, exponential amount of computational resources is 
needed, typically oc exp(ht) where h is the Kolmogorov's dynamical entropy related 
to exponential sensitivity to initial conditions. But on the other hand, how difficult is it 
to simulate isolated and bounded quantum systems of many interacting particles using 
classical resources? In analogy with the classical (chaotic) case, we might expect that 
the best classical simulation of typical quantum systems (in TL) is exponentially hard, 
i.e. the amount of computing resources is expected to grow exponentially with time. 

Even though there is no exponential sensitivity to initial conditions in quantum 
mechanics, there is a tensor-product structure of the many-body quantum state space 
which makes its dimension to scale exponentially with the number of particles, as 
opposed to linear scaling in the classical case, and due to presence of entanglement 
generic quantum time evolution cannot be reduced to (efficient) classical computation 
in terms of non-entangled (classical like) states. Here we propose the idea to use 
the computation complexity of best possible classical simulation of quantum dynamics, 
as a measure of quantum algorithmic complexity. This section essentially reviews 
the article [50J. We note that our proposal is different from existing proposals of 
quantum algorithmic complexity [511 [52J, [53J, [51] , namely we consider merely the classical 
complexity of (best) classical simulations of quantum states. In the sense of Mora and 
Briegel [54J, quantum algorithmic complexity per unit time of initially simple time- 
evolving quantum states propagated by locally interacting Hamiltonian H is clearly 
vanishing, since an approximate quantum circuit which reproduces the state after time 
t is a simple repetition of Trotter-Suzuki decomposition of exp(— iHSt). 

5.1. Time dependent density matrix renormalization group: How far can it go? 

Recently, a family of numerical methods for the simulation of interacting many-body 
systems has been developed [8] which is usually referred to as time-dependent density- 
matrix-renormalization group (tDMRG), and which has been shown to provide an 
efficient classical simulation of some interacting quantum systems. Of course, it cannot 
be proven that tDMRG provides the best classical simulation of quantum systems, but 
it seems that it is by now the best method available. Simulations of locally interacting 
one-dimensional quantum lattices were actually shown rigorously to be efficient in the 
number L of particles [55] (i.e., computation time and memory resources scale as 
polynomial functions of L at fixed t, whereas here we are interested in the scaling 
of computation time and memory with physical time t (in TL, L = oo), referred to as 
time efficiency. 

In this section we address the question of time efficiency implementing tDMRG 
for a family of Ising Hamiltonian (J7j) which undergoes a transition from integrable 
(transverse Ising) to non-integrable quantum chaotic regime as the magnetic field is 
varied. We mainly consider time evolution in operator spaces [10J, say of density matrices 
of quantum states, or (quasi) local operator algebras. Note that time evolution of pure 
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states is often ill defined in TL [51]. As a quantitative measure of time efficiency we 
define and compute the minimal dimension D e (t) of matrix product operator (MPO) 
representation of quantum states/observables which describes time evolution up to time 
t within fidelity 1 — 0(e). Our main question concerns possible scaling of D e (t) for 
different types of dynamics, and indeed we shall demonstrate a correspondence between, 
respectively, quantum chaos or integrability, and exponential or polynomial growth of 
D e (t). 

The key idea of operator-valued tDMRG pH] is to represent any operator in a 
matrix product form, 



o MP o = J2 tT ( A o A T---A 



L-l J a °\ 



a 



SL-l 

L-l > 



(16) 



in terms of 4L matrices A s - J of fixed dimension D. The number of parameters in the MPO 
representation is ALD 2 and for sufficiently large D it can describe any operator. In fact, 
the minimal D required equals to the maximal rank of the reduced super-density-matrix 
over all bi-partitions of the chain. The advantage of MPO representation lies in the fact 
that doing an elementary local one or two qubit unitary transformation O' = U^OU can 
be done locally, affecting only a pair of neighboring matrices . 

We will illustrate evolution of Ising chain (JTj), with open boundary conditions, for 
two different magnetic field values: (i) an integrable (regular) case H R = #'(1,0,2) 
with transverse magnetic field, and (ii) non-integrable (quantum chaotic) case Hq = 
H'(l, 1, 1) with45° tilted magnetic field. To confirm that Hq, and _£/r, indeed represent 
generic quantum chaotic, and regular, system, respectively, we calculated level spacing 
distribution (LSD) of their spectra (shown in fig. [3]). LSD is a standard indicator of 
quantum chaos [2]. It displays characteristic level repulsion for strongly non-integrable 
quantum systems, whereas for integrable systems there is no repulsion due to existence 
of conservation laws and quantum numbers. Evolution by tDMRG proceeds as described 




Figure 3. Nearest neighbor LSD for He (left) and (right) for L = 12. Dashed 
curves are p(s) — sn/2 exp (— 7r 2 s 2 /4) (left) and p(s) = exp(— s) (right), typical for 
chaotic and regular systems, respectively [2] . Eigenenergies € [—9,9] were used and 
statistics for even and odd parity states were combined. 



in [8], [TOl [50] using an approximate Trotter-Suziki factorization, for some time step 5t, 
of the evolution operator U(t) = exp(— iHt) in terms of 2-qubit gates. And for each 
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two qubit gates, the matrices A a ? can be updates using a singular value decomposition 
with some truncation error rj. Interestingly, it has been found [50] that the sum of all 
truncation errors up to time t, denoted by r](t) (provided St is small enough so that the 
Trotter-Suzuki factorization error is negligible) is simply proportional to fidelity error, 
namely 

l-F(t)^c V (t)/St, (17) 

where 

. \tr{0 M po{t)0 

exact 

U |tr{0^ PO (t)}||tr{O e 2 xact (t)}|' 1 > 

and numerical constant c ~ 1 does not dependent on either St, D or L. The central 
quantity we are going to study is D e (t) which is the minimal dimension D of matrices 
A}* in order for the total truncation error i](t) to be less than some error tolerance e, 
for evolution to time t. In numerical experiments shown we take e = 10~ 4 except for 
simulating thermal states of quenched Hamiltonians where e = 1CT 6 . 



5.2. Simulating pure states 

For a reference we start by investigating the evolution of pure states following the 
basic tDMRG algorithm [8]. We studied time efficiency of simulation of pure states in 
Schrodinger picture, for which many examples of efficient applications exist, however all 
for initial states of rather particular structure, typically corresponding to low energy 
sectors of few quasi-particle excitations or to low dimensional invariant subspaces. 
Treating other, typical states, e.g. eigenstates of unrelated Hamiltonians, linear 
combinations of highly excited states, or states chosen randomly in the many-particle 
Hilbert space, we found that, irrespectively of integrability of dynamics, tDMRG is not 
time-efficient, i.e. D £ (t) typically grows exponentially in time even in the integrable 
case of transverse field (consistently with a linear growth of entanglement entropy |57j). 
Numerical results are summarized in fig. HI 



5.3. Simulating local observables 

We continue by discussing the time efficiency of operator-valued tDMRG method using 
MPOs (EEEJ). Let us first study the case where the initial operator is a local operator 
in the center of the lattice O(0) = cr s L ^ 2 ,s G {x, y, z}. In the integrable case time 
evolution 0(t) can be computed exactly in terms of Jordan- Wigner transformation and 
Toeplitz determinants [HE], however for initial operators with infinite meieaQ like e.g for 
o~^j 2 , L — > oo, the evolution is rather complex and the effective number of terms (Pauli 
group elements) needed to span 0(t) grows exponentially in t. In spite of that, our 
numerical simulations shown in fig. [5] strongly suggest the linear growth D e (t) ~ t for 
initial operators with infinite index. Quite interestingly, for initial operators with finite 

X Index of a product operator [Sect. 2, 1st of Refs.[58]] is half the number of fermi operators in 
Jordan- Wigner transformation of O and is a conserved quantity for Hr. 
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Figure 4. MPS rank D e (t) for simulating pure states with integrable transverse Ising 
model -ffpt, except full squares which are for Heisenberg XX chain, starting from the 
initial states indicated in the legend (explanation: \W) = (|10 . . . 0) + |01 . . .0) + 
•••|00...1))/VIand \GHZ) = (|00 . . . 0) + |11 . . . 1))/V2). Note that the full squares 
corresponds to the same example as studied in Ref. [S] • 

index, D e (t) saturates to a finite value, for example D t (oc) = 4 for o- z L j 21 or D e (oo) = 16 
for 0£/2-i°l/2- m non-integrable cases the rank has been found to grow exponentially, 
D e (t) ~ exp(/iqt) with exponent h q which does not depend on e, properties of O(0) or 
L, for large L. For H = Hq we find h q = 1.10. 

5.4- Simulating extensive observables 

In physics it is often useful to consider extensive observables, for instance translational 
sums of local operators, e.g. the Hamiltonian H or the total magnetization M s = 
J2jZo a j- As opposed to local operators, extensive initial operators, interpreted as W- 
like states in operator space, contain some long-range 'entanglement' so one may expect 
that tDMRG should be somewhat less efficient than for local operators. Indeed, in 
the integrable case we find for extensive operators with finite index that D e (t) does 
no longer saturate but now grows linearly, D e (t) ~ t, whereas for extensive operators 
with infinite index the growth may be even somewhat faster, most likely quadratic 
D e (t) ~ t 2 but clearly slower than exponential. In the non-integrable case, we again 
find exponential growth D e (t) ~ exp(h q t) with the same exponent h q = 1.10 as for local 
initial observables. The results are summarized in fig. [6l 

5.5. Simulating thermal states 

In the last set of numerical experiments we consider time efficiency of the evolution of a 
thermal initial state 0(t) = Z~ l exp(—f3Ho) under a sudden change of the Hamiltonian 
at t = 0, namely H(t < 0) = H = H'(l,0,l),H(t > 0) = H v Again, we treat 
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Figure 5. D e (t) for local initial operators. We consider three cases O(0) = cr^j^ 
(empty circles, squares and triangles), for non-integrable evolution He, and four cases, 
O(0) = cr*'J 2 (full squares, diamonds), cr £,/2-i cr L/2 (^ u ^ triangles), all with infinite 
index, and O(0) = < j l/2-i' j l/2 (^ uu cn " c l es ) with index 2, for integrable evolution Hr. 
In the inset we plot the data for the non-integrable case He in semi-logarithmic scale, 
and the full line in the inset illustrates exponential growth oc e 11 *. Full squares and 
diamonds are for L = 40, otherwise L = 20. 




Figure 6. D e (t) for extensive initial operators. For both Hamiltonians He, H-r, we take 
O(0) = a* (empty, full squares) with infinite index, and O(0) = H'(l, 0, 1) (empty, 
full circles) with index 1. For H-r we also show the case O(0) = a j a j+i + ^i^J+i 
(full diamonds) with index 1 and 2. In the semi- log inset we illustrate exponential 
increase oc e 11 * (full straight line) for Hq and polynomial ~ t 2 (full curve) for Hr. 
For full circles L — 64, otherwise L — 32. 



two situations: in the first case we consider change after which the Hamiltonian remains 
integrable, H 1 = H'(\, 0, 2) = H R , while in the other case the change breaks integrability 
of the Hamiltonian, Hi = H'(\, 1, 1) = H c . The initial state is prepared by means of 
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Figure 7. D e (t) for thermal states of Ho with f3 = 0.01 {(3 = 0.05 in inset), for 
evolution with He (open symbols) and iJpt (full symbols) at L — 40. Solid curves 
again indicate exponential increase oc e 11 '. 

imaginary time tDMRG from initial identity operator using the same MPO rank D as 
it is later used for real time dynamics. We find, consistently with previous results, that 
at high temperature (|3 < 1) the rank D e (t) grows very slowly, perhaps slower than 
linear, in the integrable case, and exponentially D e (t) ~ exp(h q t), in the non-integrable 
case. Interestingly, at lower temperatures we find exponential growth in both cases, 
even in the integrable one. This is not unreasonable as the initial (thermal) state can be 
expanded in a power series in f3 and the higher orders Hq become less local with longer 
entanglement range as we increase the power p. These results are summarized in fig. [7J 

6. Quantum chaos and far from equilibrium quantum transport 

In this section we would like to demonstrate the connection between quantum chaos, 
(non)integrability and transport in non-equilibrium steady states of interacting quantum 
chains [59]. Within the linear response theory, the property of quantum mixing (which is 
typically implied by quantum chaos, i.e. by the validity of random matrix level statistics 
|24j . see also correlation functions in sect HD , is typically synonymous for normal quantum 
transport since it implies finite Kubo transport coefficients [60] (provided temporal 
correlation functions decay fast enough). 

However, here we would like to address the connection between quantum chaos and 
transport in far-from-equilibrium steady state, which may be beyond applicability of 
linear response. In particular, we are interested in the validity |61j of the Fourier law 
J = — kVT in quantum chains, relating the macroscopic heat flux J to the temperature 
gradient VT. 

To investigate this problem one has to deal with a finite open system connected to 
heat baths. Here we consider an interacting quantum spin- 1/2 chain (J7J) which exhibits 
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the transition from integrability to quantum chaos as a parameter, e.g. the magnetic 
field, is varied. The standard treatment of this problem is based on the master equation, 
thus limiting numerical investigations to relatively small system sizes. By using this 
method, in an interesting paper [62J, the decay of current correlation function in a 
model of non-integrable chain of quantum spins is computed. However, these results 
were not fully conclusive and the conclusions rely on linear response theory. Also, in 
[63J Lindblad formalism was used to study the validity of Fourier's law for different type 
of spin-spin interaction. 

Here we describe a different approach (see Ref.[59] for details) namely we follow 
the evolution of the system described by a pure state, which is stochastically coupled 
to an idealized model of heat baths. Stochastic coupling is realized in terms of a local 
measurement at the boundary of the system and stochastic but unitary exchange of 
energy between the system and the bath. By this method we have been able to perform 
very effective numerical simulations which allow to observe a clear energy/temperature 
profile and to measure the heat current J. Again we consider Ising spin chain in the 
magnetic field ([7]) of size L, where the first and the last spin are coupled to thermal 
baths at temperatures Tj and T r , respectively. In the non-integrable regime where the 
spectral statistics is described by RMT - the regime of quantum chaos - we found very 
accurate Fourier law scaling J /AT oc 1/L, where L is the size of the chain. In the 
integrable and near-integrable (non-ergodic) regimes instead, we found that the heat 
transport is ballistic J oc L°. 

Let us describe the numerical simulations. Again we consider the autonomous 
model ([7]) at three particular cases: (i) quantum chaotic case H = H'(—2, 2, 3.375) at 
which LSD agrees with RMT and thus corresponds to the regime of quantum chaos, (m) 
integrable case H = H'(—2, 0, 3.375), at which LSD is close to the Poisson distribution, 
and (Hi) intermediate case H = H'(—2, 2, 7.875) at which the distribution LSD shows 
and intermediate character of weak level repulsion and exponential tail [59J. 

Let us now turn to study the energy transport in this model system. To this end we 
need to couple both ends of the chain of spins to thermal baths at temperatures T T ,T\. 
We have devised a simple way to simulate this coupling, namely the state of the spin 
in contact with the bath is statistically determined by a Boltzmann distribution with 
parameter T. Our model for the baths is analogous to the stochastic thermal baths used 
in classical simulations [M] and we thus call it a quantum stochastic thermal bath. In 
the representation basis of the wave function at time t can be written as 

m))= E C S0!S1 _ SL _ 1 (t)Us 1 ...s L _ 1 ) , (19) 

s ,si,...,sl-i 

where s n — 0, 1 represents the up, down state of the n-th spin, respectively. The wave 
function at time t is obtained from the unitary evolution operator U(t) = exp(-iHt). 
The interaction with the bath is not included in the unitary evolution. Instead, we 
assume that the spin chain and the bath interact only at discrete times with period r at 
which the states of the leftmost (sq) and the rightmost (sl-i) spins are stochastically 
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reset. Thus, the evolution of the wave function from time t to time t + r can be 
represented as 

mt + T)) = ~(f3 h (3 r )U(T)\4>(t)) , (20) 

where S(/?i, /3 r ) represents the stochastic action of the interaction with the left and right 
baths at temperatures /3j and (3~ l respectively. 

The action of E(/3\, /3 T ) takes place in several steps: 

• (i) The wave function is first rotated by the angle a = arctan(/i x //i z ) to the 
eigenbasis of the components a\ = h ■ Bo/h, o x = h ■ Bl-x/H of the edge spins 
along the field h = (h x , 0, h z ), that is \ip) -> e - ia(CT « +,T ^i )/2 |^). Here, /i = \h\ stands 
for the magnitude of the magnetic field. 

• (m) A local measurement of the observables o~i, o r is performed. Then the state of 
the spins at the borders collapses to a state (so> s !-i) with probability 

pO»S,«U)= E l^-s..!,...,.^-! I 2 ■ ( 21 ) 

So, after choosing Sq,s* l _ 1} we put all coefficients C SOjSlv .. iSi l with (so,sz-i) 7^ 
(sq, s£-i) to zero. 

• (Hi) The new state of the edge spins (s ,s,l_i) is stochastically chosen. After this 
action simulating the thermal interaction with the baths each of the edge spins is 
set to down, (up) state with probability — /x). If the new state of sq, or sl-i, 
is different than the one after step (ii), then a simple unitary spin flip is performed 
to the wave-function. The probability fi(/3) depends on the canonical temperature 
of each of the thermal baths: 

= e-fif+efiH ; i^ftr}. (22) 

• (iv) Finally, the wave function is rotated back to the o~ z n basis, \ip) — > 

This completes the description of the interaction with the quantum stochastic bath. 
This interaction thus (periodically) resets the value of the local energy ho~i jT of the spins 
in contact with the baths. Therefore, the value of r controls the strength of the coupling 
to the bath. We have found that, in our units, r = 1 provides an optimal choice. We 
have nevertheless performed simulations for other values of r with qualitatively similar 
results. In particular, for weak couplings (r > 1) the heat conductivity does not depend 
on the coupling strength . 

Note that our method does not correspond to a standard stochastic unraveling of 
a master equation for the density operator (e.g., in the Lindblad form) [65]. However, 
we have tested that, using our prescription, averages over the ensemble of quantum 
trajectories or time averages of one given quantum trajectory are sufficient to reconstruct 
a density matrix operator p = \tp(t))(tp(t)\ that correctly describes the internal thermal 
state of the system in and out of equilibrium. For each run the initial wave-function 
|^(0)) of the system is chosen at random. The system is then evolved for some 
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Figure 8. Out of equilibrium energy profile (H n ) for the chaotic chain. The 
temperatures of the baths T\ = 5 and T r = 50, are both in the high temperature 
regime. Results for chains of size L = 8 (crosses), L = 14 (open circles) and L = 20 
(solid circles) are shown. The dashed line was obtained from a linear fit of the data 
for L = 20 for the L — 4 central spins. Insets (I) and (II) show the energy profile for 
the integrable and intermediate cases respectively, for L = 15. 



relaxation time t ic \ after which it is assumed to fluctuate around a unique steady state. 
Measurements are then performed as time averages of the expectation values of suitable 
observables. We further average these quantities over different random realizations of 
"quantum trajectories" . In order to compute the energy profile we write the Hamiltonian 
<^ in terms of local energy density operators H n : 

H n = Ja*a* +1 + - ■ (a n + a n+1 ) . (23) 

such that the total Hamiltonian ([7]) can be written as H = ^ n H n apart from the 
boundary corrections. 

First we have performed equilibrium simulations in order to show that time averaged 
expectation values of the local energy density can be used to determine canonical local 
temperature. To this end we set the left and right baths to the same temperature T. 
For low T, the energy per site E = (1/L)(H) saturates to a constant which, together 
with the entire energy profile (H n ), is determined by the ground state. However, for 
larger T > 1, the energy profile is constant within numerical accuracy, and numerical 
simulations give E ~ —1/T, all results being almost independent of L for L > 6. The 
numerical data for E(T) can be well approximated with a simple calculation of energy 
density for a two-spin chain (L = 2) in a canonical state at temperature T, namely 
E ca , n (T) = trHoe~ H °/ T /tre~ H °/ T . Therefore, if the temperatures of both baths are in 
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Figure 9. Size dependence of the energy current in the chaotic chain with T\ — 5 and 
T T = 50. We show J/ AE (open circles) and J/ AT (open squares). The dashed lines 
corresponds to 1/AL scaling. In the inset, the size dependence of the energy current 
is shown for the integrable (solid circles) and the intermediate (solid squares) cases. 



high T regime, then we can define the local temperature via the relation T oc — l/E. We 
stress that equilibrium numerical data shown are insensitive to the nature of dynamics 
(consistent with results of Ref.|66j), whether being chaotic, regular or intermediate. 

In fig. [8] we show the energy profile (H n ) for an out of equilibrium simulation of the 
chaotic chain. In all non-equilibrium simulations, the temperatures of the baths were 
set to T\ — 5 and T r = 50. After an appropriate scaling the profiles for different sizes 
L collapse to the same curve. More interesting, in the bulk of the chain the energy 
profile is in very good approximation linear. In contrast, we show that in the case of 
the integrable (inset I) and intermediate (inset II) chains, no energy gradient is created 
which is a characteristic of ballistic transport. 

We now define the local current operators through the equation of continuity: 
d t H n = i[H, H n ] = -( J n+1 - J n ), requiring that J n = -i[H n , # n _i]. Using eqs. ([23]) and 
([7]) the local heat current operators are explicitly given by J n = h z J (cr*^ — cr^ +1 ) o y n . 
In fig. [9] we plot J / AE as a function of the size L of the system for sizes up to L = 20. 
The mean current J is calculated as an average of (J n ) over time and space n. The 
energy difference was obtained from the energy profile as AE = (Hl^ s ) — (H 2 ). Two 
spins near each bath have been discarded in order to be in the bulk regime. Since 
AL = L — 5 is an effective size of the truncated system, the observed 1/AL dependence 
confirms that the transport is normal (diffusive). Moreover, also the quantity J/ AT, 
where AT = —1/(Hl_ s ) + 1/(H 2 ), shows the correct scaling with the size L. On the 
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other hand, in integrable and intermediate chains we have observed that the average heat 
current does not depend on the size J oc L°, clearly indicating the ballistic transport. 

7. Quantum relaxation and complexity in a toy model: Kicked Ising Chain 

Let us now come back to the kicked Ising model (JTJ) and try to consider some very 
elementary but fundamental questions considering its dynamics and non-equilibrium 
statistical mechanics. Observing data of fig. [I] in section (j3J) one can conclude that the 
model perhaps displays an interesting order to chaos, or non-ergodicity to ergodicity 
transition when the integrability breaking parameter is increased. Now we would like to 
inspect this transition more closely, and in particular understand the rate of relaxation 
to equilibrium in the ergodic and mixing case. We should stress right from the start 
that we are unable to prove any non-trivial statements about the model, but we can 
provide many suggestive numerical experiments which can be performed in a rather 
efficient way. We have learned in section [5] that it may be more fruitful to consider time 
evolution in the operator algebra spaces instead of in the spaces of pure states. Let us 
go now a bit deeper into this subject. 

For some related results on high-temperature relaxation in isolated conservative 
many-body quantum systems see e.g. Refs.[67J EH EH [TQl [7T] . 

7.1. Time automorphism 

Time automorphism on unital quasi-local C* algebra (see e.g. [72] for introduction into 
the subject) Az, T : Az — > Az of an infinite KI lattice, for one period of the kick, can 
be explicitly constructed by the following observations. 

Formally, for any A e Az, TA := WAU, where U is given by either (!3f4f6|) . Let 
A[ mtn ], with m < n, denote a finite, local 4 n - m+1 dimensional algebra on a sub-lattice 
[m,n] C Z, which is spanned by operators G^a^+i • • It is straightforward to 

prove that dynamics is strictly local 

T : A[ m ,n] — ► A[ m -l,n+l]- (24) 

In other words, the homomorphism T (124p is a simple nontrivial example of a quantum 
cellular automaton as defined by Schumacher and Werner [56] . 

Az can also be treated as a Hilbert space with respect to the following inner product 
(A\B) = u(A*B), where u(A) is a tracial state, uo(A) = 2-^- m+1 hrA for A e A [m , n] . 
This Hilbert space can be in fact considered as a ID lattice of 4- level quantum systems 
(qudits with d = 4) with the orthonormal basis | . . . , s_i, s , s 1; ...) = •• • ct^Vq ^ 1 • • • 
labeled by an infinite sequence of 4-digits Sj G Z 4 . Restricting for a moment to a dimer 
lattice A[jj + i] we can write the adjoint action of a 2-qubit gate W (JSJ) in terms of a 
16 x 16 unitary matrix 

W j,j+i\sj,s J+ i) = W^a/a^W = l r ^ r i+i) w (r 3 ,r 3+1 ),(s 3 ,s J+1 ) (25) 

rj,r j+1 eZ 4 
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where, very explicitly 

W ( , 1)r2)>(si)S2) = ±tr [K 1 ® o*») WV* ® a S2 )iy] . (26) 

We should note that the map is unital W|0, 0) = 1 0, 0) , and that due to anti-unitary 
symmetry of KI model, the matrix W is real. We extend the map Wj J+1 to entire 
algebra A% by W j:j+1 (A ® B) — Wjj+^A) <g> B, for any A G Ayj+i], B G Az~\j,j+i]- 

Now, following the protocol ([6]) we finish the construction of the time automorphism 
as a string of right-to-left ordered 2-qudit (with d = 4) gates 

Explicitly, for any local observable A = E Sm , Sm+1 ,..., Sn a( Sm , Sm+1 ,..., Sn )|s m , s m +i, ■ ■ ■ , O e 
-4[m,n]) we have the following 
Algorithm 1: 

(i) Set an initial vector: a(°i_ 1)Sm! ..., Sn , Sn+l) = <S« m -i,ofl(» rai ...,,«)<S» n+ i,o- 

(ii) For = 0, 1, . . . , n — m + 1: 

^= V W (s , its .^frwiaf , v (28) 

(s m _i,s m ,..,s„,s„ + i) / j \s m -l+k>Sm+k)A r > T ) {s m - 1 ,...,s m+k _ 2 ,r,r',s rn+k+1 ,...,s n+1 ) \ > 

r,r'gZ4 

(hi) The result is TA = ^"Tf s , «n+i)- 

The algorithm produces exact result in (n — m + 2)4" _m+4 multiplications and about 
the same number of additions. Let P[ m , n ] : Az — > A[ m , n ] denote a linear orthogonal 
projector, satisfying P^ m>n ^(A <g> E) = (t\E)A if A G ^4[ miTl ],£? G Az-[ m ,n]- Let us define 
a truncated time evolution operator T[ m>n ] = P[ TO>n ]T : A[ m ^ n ] — > *4.[ m ,n] which we can 
actually implement on a computer with a finite memory register. 

The "infinite-temperature" time correlation function of (traceless) local quantum 
observables A,Be A[ m , n ] can be written as Cba^) = {B\T t A), where f £ Z. An 
interesting question is, up to what time t the CBAif) can be computed numerically 
exactly with a finite computer register [m — I, n + 1] of r = n — m + 1 + 21 qudits of size 
4 r ? Due to locality (1241) of time homomorphism one can easily prove that 

(B\T t A) = (B\T\ m _ lin+l] A), for t < 21 (29) 

hence the correlation functions are computable exactly up to time 21, and as we shall see 
later, truncated correlation function C^' A (t) = (B\T^ m _ rn+r ^A) often well approximates 
CsAit) even at later times, or even its asymptotic decay. 

Let us continue our discussion by considering time evolution for translationally 
invariant extensive (TIE) observables. Given some quasi-local observable A G Az 
we shall construct the corresponding TIE observable by a formal mapping, A — > 
F(^4) = XLez^z^) in terms of lattice translation automorphisms S^, : Az — ► Az, 
S x ((Tj) = o~j +x . The image of the entire quasi-local algebra under this mapping 
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Z = F(Az), is not a C* algebra, but it is a linear space which can be again turned 
into a Hilbert space with the following inner product 

{{X\Y)) = Jim -^(p { _ n>n] ( X )\P hn>n] (Y)) (30) 

where the domain of projector Pr ntTl \ is extended to Z by continuity. Orthonormal 
basis of Z is given by TIE observables ^( Co ,ci,...,c r _i) = F^Vi 1 • • • o~ c r r Si), for orders r = 
1,2,.. ., and for uniqueness of notation, requiring c , c r _i ^ 0. We shall interchangeably 
represent finite sequences of 4-digits with integers, (c , Ci, . . . c r _i) = c = ^j=o c j^ ■ Let 
Z r be 3 x 4 r_1 dimensional subspace spanned by TIE observables Z c with order < r, i.e. 
for c having at most r base-4 digits, so we have an inclusion sequence Z\ C Z 2 . . . C Z . 

Since time and space automorphisms commute TS X . = S X T, one can immediately 
extend the time map onto the space of TIE observables, T : Z — > Z by continuity. 
Formally, we have TF = FT. Furthermore, locality (J24l) implies 

f : Z r -> Z r+2 , (31) 

so we can write a simple adaptation of Algorithm 1 for explicit construction of a time 
map of an arbitrary finite-order TIE observable Y = 2^o<c<™° d ^ Uc^c £ Z r : 
Algorithm 2: 

(i) Take the following pre-image of the TIE observable A = Y,t< c <T A ^ ^l 4c ) e Ahr\, 
namely F(A) = Y. 

(ii) Compute a c of T(A) = J2o< c <4 r + 2 °cl c ) ^ ^[o,r+i] according to Algorithm 1. 

(iii) Transforming back to Z the result reads 

c^o (mod 4) a c + a Ac + ai6c, c < 4 r ; 

f (Y) = F(T(A)) = y'c Z " y'c={ *c + a Aci ¥<c<A r+1 ; (32) 

Let us further define the natural truncations of TIE space to order r, P r : Z — > 2 as 
orthogonal projections P r (X) = 2^o<c<™° d Zc((Z c \X)), and truncated time evolution 
operators T r = P r T : 2 r — > 2 r , which are naturally implemented on a computer by 
simply truncating overflowing coefficients y' c . 

Physically interesting question now concerns computation of time correlation 
functions between a pair of finite order (say q) TIE observables X, Y e Z q , namely 
Cyx(t) = ((y|T'X)), for example in fig. [1] we have shown the case of X = Y = Z%. As 
a consequence of locality (|24|) . and translational invariance, we find that the truncated 
evolution on Z r , reproduces correlation functions exactly 

((BIT' A)) = p|T* A)), for t < r - q. (33) 

Few remarks are in order: (i) Truncated translationally invariant time evolution T r is 
perhaps more natural object to study than truncated local time evolution Tr mn i, for 
a simple reason that the truncation P r commutes with a shift S x , while P[ mj „] does 
not. (ii) A space Z can be identified with translationally invariant linear functionals 
over A z , namely X(A) = ((X\FAj), X e Z, A e A z . We have (f^X)(A) = X(TA) 
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and (TX)(A) = X(T^A) where Hermitian adjoint maps Tt and simply correspond 
to time reversed dynamics, (iii) Convex subspace of positive translationally invariant 
functionals W C Z is an interesting invariant subspace of physical states, TW C W. 

7.2. Relaxation and quantum Ruelle resonances 

In classical mechanics of chaotic systems one typically observes that states (phase- 
space densities) develop small details in the course of time evolution at an exponential 
average rate. Consequently, introducing a small stochastic noise of strength e to Perron- 
Frobenius operator (PFO, i.e. Liouvillian propagator for discrete time dynamics), makes 
it non-unitary and shifts its spectrum inside the unit circle. Typically, the effect of noise 
is equivalent to an ultraviolet cutoff - truncation of PFO - at the Fourier scale k ~ 1/e, 
and often the leading eigenvalues of truncated PFO - the so-called Ruelle resonances - 
remain frozen inside the unit circle in the limit e —>■ (or k —>■ oo) [73J. For a general 
introduction to relaxation phenomena in classical Hamiltonian dynamics see e.g. [73] . 

Let us now draw some some analogies with our quantum setting. We have seen 
that that the evolution T somehow most closely resembles Liouvillian evolution of 
classical Hamiltonian dynamics. In Hilbert space topology, operator T is unitary and its 
spectrum lies on a unit circle, just like in the case of classical PFO. However, truncated 
(3 x 4 r_1 ) x (3 x 4 r_1 ) matrices T r represent natural "ultraviolet" cutoff truncations 
for increasing orders r. Let us check numerically if some eigenvalues of these matrices 
remain frozen when r — ► oo. 

Indeed, as we demonstrate in fig. [101 we find several eigenvalues which converge 
as r increases in the case of strongly non-integrable (quantum chaotic) case, with a 
gap between an eigenvalue of maximal modulus and the unit circle, whereas in the 
integrable case a set of r eigenvalues touch the unit circle (actually eigenvalue 1 is r- 
fold degenerate). Numerical results suggest the following speculative conclusions. Let 
e~ Qn be the converged (frozen) eigenvalues of T r , and {©^}, {©^} the corresponding 
right and left eigenvectors, respectively. Then for arbitrary pair 1,7 6 2, the time 
correlation function can be expressed in terms of spectral decomposition (see e.g. [74"] ) 

r m „, (WM^IX)) 

C YX (t)~^ Wn e , w n = . (34) 

The above relation is the contribution of the point spectrum and is exact if the spectrum 
is pure-point. However, in classical cases one may quite typically have various singular 
components and branch cuts [71]. Note that the denominator ((0^|@^)) is finite although 
both vectors should have infinite I 2 norm ((0£|G£)) = oo, ((©^|©^)) = oo, for any 
eigenvalue away from the unit circle, Re q n ^ 0. 

There is a simple relation between the spectrum of PFO and the ergodic properties 
of dynamics: (i) If there is a spectral gap, i.e. there exists A > such that for all n, 
\e~ q,l \ < exp(— A) < 1, then dynamics is exponentially mixing, Cyx(t) < C exp(— Xt). (ii) 
If some eigenvalues are on the unit circle, meaning that the corresponding eigenvector 
coefficients should be in I 2 , then the system is non-mixing since there are correlation 
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Figure 10. The spectra of truncated transfer operators TV, for r = 5, 6,7 in 
strongly non-integrable case J = 0.7, /i x = 0.5, ,h z — 1.1 (left) and integrable case 
J = 0.7, h x — 0.0, h z = 1.1 (right), lying inside complex unit circle (thin arcs). 
The points in upper/lower unit semi-disks correspond to positive/negative parity 
^(co,ci,...,ci„i) — ^(ci_i ci,c ) eigenvectors. Arrows point at converged positions 
of the leading eigenvalue e~ qi . 



functions which do not decay, (iii) If some eigenvalues are at 1 then the system is 
non-ergodic since the correlation functions may have non- vanishing time-averages. If 
Q n is a complete set of orthonormalized eigenvectors corresponding to eigenvalue 1, 
(iQnlQmj) — fin,m (and note that since we are on the unit circle: = Q^) then 

D x :=C^t) = Y,\iX\Qn))\ 2 . (35) 

n 

The latter (iii) happens in generic completely integrable quantum lattices, where Q n 
correspond to an infinite sequence of conservation laws [75]. Furthermore, we have 
a strong numerical evidence that also in certain non-integrable quantum lattices [76J, 
and also in KI model [17J, one has a regime where few normalizable ('pseudo-local') 
but not local (like in integrable models) conservation laws exist. This situation we call 
the regime of intermediate dynamics and is characterized by a non-vanishing stiffness 
Dx 7^ signalling ballistic transport. 

In figure [IT] we compare the time autocorrelation function of the transverse 
magnetization M = Y^jez a j = ^3' computed in three different ways: (1) from 



Chaos and Complexity of quantum motion 



28 




10 20 30 40 50 60 70 

t 

Figure 11. Correlation function of the transverse magnetization C(t) = ({M\M(t))), 
in the mixing case J — 0.7, h x — 0.5, h z — 1.1, computed from finite system dynamics 
for different sizes L (symbols), and from truncated adjoint propagators T r of infinite 
systems (curves) for different truncation orders r. The chain line indicates the 
asymptotics based on the leading quantum Ruelle resonance. 



exact time evolution C^it) = ^(MU^MU^) on a finite lattice of length L with 
periodic boundary conditions, (2) iteration of truncated TA matrix on infinite lattice 
C r (t) = ((M|T*M)), and (3) asymptotics based on (few) leading eigenvalue resonance(s) 
(using formula fl34"l) in terms of q n and w n .). 

We note that the leading eigenvalue and eigenvector of truncated TA T r can most 
efficiently be computed using our Algorithm 2 as a key step of an iterative power- 
method. In this way we were able to perform calculations of the leading Ruelle resonances 
up to r = 15 in contrast to full diagonalization of truncated matrices T r which were 
feasible only up to r = 7. Let us observe the structure of the eigenvector coefficients 
v c' R = ((^c|©^' R )) corresponding to the leading eigenvalue. Numerical results (see fig.[T2| 
see also subsect EH later) strongly suggest self-similar behaviour upon multiplying the 
code c by 4 which is a consequence of the fractal structure of the transfer matrix (see 
illustration in Ref . [771 [18] ) . 

The stiffness Dx and the spectral gap A = 1 1 — e~ qi | may be considered as order 
parameters, characterizing a particular kind of dynamical phase transition, namely the 
transition from non-ergodic dynamics - ordered phase, where Dx ^ for a typical A0, 
or A = 0, to an ergodic and mixing dynamics - disordered phase, where Dx = 0, for 
all traceless X e Z, or A > 0. We know that KI model is non-ergodic in integrable 
regimes. Let us consider a fixed transverse field case h x — 0, and start to switch on a 
small amount of longitudinal field h x . The interesting question is whether the transition 



§ Any observable X which is not orthogonal to all conservation laws Q ni i.e. normalizable eigenvectors 
of T with eigenvalue 1. 
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Figure 12. Eigenvector Qf- — J2 c Vc ^ c ^ °f the leading eigenvalue (closest to unit 
circle) has statistically selfsimilar structure, when expanded in Z c . We plot the 
modulus of coefficients v c (in log scale) of the right eigenvector versus the log (in base 4) 
of the integer code c. Dashed line indicates power law scaling c~ v with slope v = 0.32. 
In the inset we plot partial scalar products u m — ^^°-i^°^ i t!_ 1 ({G>\\Z c ))((Z c \Q^')) 
with the corresponding left eigenvector within fixed orders m. 
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Figure 13. Two dimensional numerical phase diagram for kicked Ising lattice at cutoff 
order r = 7. Gray level indicates the spectral gap log A of T r as a function of h x and 
h z at fixed J = 0.7. White regions correspond to A < 10~ 6 . 
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Figure 14. Type I transition: Spectral gap of T r as a function of h x at J = 0.7, h z — 
1.1 at different truncation orders r. In the inset we indicate a line of transition in 2D 
phase diagram. 
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Figure 15. Type II transition: Type I transition: Spectral gap of T r as a function of 
/i x = h z at J = 0.7 at different truncation orders r. In the inset we indicate a line of 
transition in 2D phase diagram. 
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happens for infinitesimal integrability breaking parameter h x in TL, or at a finite - 
critical field. In fig. [13], we fix J = 0.7 and plot a two dimensional phase diagram of 
the spectral gap A(h x ,h z ). It is clear that we have different behaviours in different 
regions of parameter space, for example we identify two: (i) Type I transition: if the 
transverse field is roughly on the interval h z G [0.7,1.2], then the spectral gap opens 
in the fastest possible manner which is allowed by a h x — > — h x symmetry and the 
analyticity of the problem, namely A oc h x . See fig. [HI (ii) Type II transition: if the 
initial transverse field h z < 0.7, or h z > 1.2, then the gap opens up in a much more 
abrupt - perhaps a discontinuous way. We give an example by scanning the diagonal 
transition, i.e. putting h x = h z and increasing h x from zero. Numerical results, shown in 
figfl5l cannot be made fully conclusive, but they are not inconsistent with a conclusion 
that an abrupt transition to ergodic behaviour takes place at h x = h z « 0.3. 

7. 3. Translationally invariant conservation laws as matrix product operators 

There is another possibly interesting way of characterizing the transition, i.e. in terms 
of pseudo-local translationally invariant conservation laws [76]. Such conservation laws 
are the square normalizable elements Q G Z, which are mapped onto themselves under 
the dynamics T(Q) = Q. We shall first make a non-trivial variational MPO ansatz for 
elements of Z, namely let us take an auxiliary vector space C D = C Dl © C° 2 © C° 3 , 
where D = Di + D 2 + D 3 . Then any operator Q, which is formally written in terms of 
MPO on the infinite spin chain 

Q= (ZL---A s -iA S0 A s ^--a R )---a s _- 1 1 a s °a?--- (36) 

...s_isosi...GZ4 

where A s G C DxD , a^, cir G C d have the block matrix form (for k = 1, 2, 3): 

/l \ /0 * * \ /*\ /0\ 

A°= 1 , A k = , a L = , a R = * , (37) 

\o o eJ \o * eJ \oJ \oJ 

represents a translationally invariant pseudo-local operator, i.e. an element of Z, 
provided that ||-Eq|| < 1 and ||-Ei|| 2 + H-E2II 2 + II-S3II 3 < !> where ||.|| is a spectral 
matrix norm and *'s stand for arbitrary matrices/vectors. Of course, converse cannot 
be generally true, not any element of Z can be written as MPO (|36|) with finite D, but 
still there are elements of the form (1361 1371) which are not in Z r for any finite r. 

There exist a straightforward algorithm which performs time evolution on MPO 
data fl37|) . namely T(Q) is also of the form (I36II37P with dimension D' < 2D. We shall 
now make the following simple numerical experiment. Let us fix D, setting D\ = D2 = 1 
representing the simplest T-invariant subclass of (I36ll37p . and optimize (maximize) the 
fidelity-like quantity 

Fln , _ l|g|Tg|| 

FiQ) - "loioT' (38) 

within this class of operators. Let us write an operator which maximizes F(Q) for a 
given D as Qe>. Note that 1 — F(Q D ) gives a strong-topology measure of conservation of 
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Figure 16. Optimized fidelity (|38jl - in log scale - for approximate conservation laws 
within MPO spaces of different fixed dimensions D (indicated in the legend) for an 
infinite KI lattice along the diagonal type 11 transition with h x = h z and J = 0.7. A 
simple stochastic search has been used to maximize fidelity F(Q). Note the similarity 
with the gap curve shown in fig ll5l 

observable Qd in one step of time evolution, so F(Qd) = 1 only for exact conservation 
laws. Increasing D may improve fidelity, if pseudo local conservation laws exist to which 
Qd may converge, however in ergodic and mixing situation where no exact pseudo-local 
conservation laws exist, increasing D should have no significant effect to fidelity F(Qe>). 
This is exactly what we observe in KI model following a line of type II transition (see 

fig. USD. 

7.4- Operator-space entanglement measures and complexity of time evolution 

Numerical results of section [5] suggested that operator space entanglement measures can 
be used to characterize the complexity of time evolution, namely the minimal required 
rank D € of MPO ansatz is simply related to entanglement entropy of a time-evolving 
local observable, which is interpreted as a Hilbert space vector. In order to make things 
as simple and precise as possible we go back to the time evolution automorphism T 
over the quasi-local spin algebra A%. Let us take some truncation order r = In + 1 
and consider the truncated map T Starting with local operators on a single site 
-4.[o,o] this truncated map is exact up to time t = n. Writing time evolved observable 
at any instant of time as A(t) = T|_ nra jA(0) = a[*_ n s Js_n, . . . , s n ) in terms of 
a 'wave-function' s and partitioning a sub-lattice at m, —n < m < n, as 

[— n, n] = [—n,m — 1] U [m, n], we can define 4 n ^ m+1 x 4 n - m + 1 reduced super-density 
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Figure 17. Purity p(°' n )(i) of reduced operator space density matrix for KI lattice, 
with different cutoff sizes r = 2n + 1 = 1 5, 1 3, 11, and for two initial centered local 
operators X = ctq and Z = erg. Note exponential decay of purity for QC case and 
saturation or slow decay of purity for non-ergodic cases (IN, NE). 

matrix as 

(s m ,,...Sn)}(Sy rt )-"? s n) ' J \S—rif--,Sm—l , s mv s n) (S— n v ■ - ; s m — 1 > s m v 5 ™) 

Since A(t) is interpreted as a 'pwre state', namely a vector from A[- n , m -\] ® A[ m , n ], 
the operator space entanglement is most simply characterized either by Von Neuman 
entropy S^ n \t) = -tiR^ n HogR^ m ^ or linear entropy S^ m ' n \t) = - log P( m ' n \t) 
where p( m ' n )(t) = tr[i?( m ' n ^(t)] 2 is a purity of reduced super-density matrix. 

In fig. [T7] we plot the entanglement purity P^°' n \t) - for close to symmetric 
bipartition m = where the entanglement is expected to be maximal - for three different 
characteristic cases of KI model, which will be in the following referred to as: quantum 
chaotic (QC), J = 0.7, h x = 0.9, h z = 0.9, mtegrable (IN), J = 0.7, h x = 0, h z = 0.9, 
and non-ergodic (NE) non-integrable case, J = 0.7, h x = 0.2, h z = 0.2. We find that 
in QC case purity decreases exponentially p(°' n \t) = exp(—h q t), meaning S^t) = h q t, 
where the exponent h q is independent of the initial observable A(0) and asymptotically 
independent of r. On the other hand, in IN case P(°' n \t) does not decay at all so the 
resulting dynamical entropy h q = 0, whereas in NE case p(°' n \t) decays slowly, likely 
slower than exponentially. 

7.5. Scaling invariance and the problem on semi-infinite lattice 

The eigenvectors of T corresponding non-unimodular eigenvalues seem to exhibit a 
certain scaling invariance (fig. [T2l . Here we would like to explore this property in 
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a little bit more detail. For that purpose we again explore the map Ti^j on the 
quasi-local algebra A]- nn ] since the representation of dynamics is conceptually simpler 
(Algorithm 1) than dynamics T2 n +i on -Z^n+i- 111 particular, it is worth to mention that 
if one traces out an additional qudit after each time step, then the dynamics is exact on 
a closed set of 4 n - m+1 x 4™~ m + 1 super-density matrices and has a simple explicit form 
in terms of a completely positive matrix map [Jj] 

R ( m+ l,n+l) {t +1)= tr(j | Tn _ m ^M ({) g, S()o ] Tj t_ m | ; (4Q) 

T r = n if ®^®if r - fc) , (4i) 

0<fc<r 

namely no truncation is needed since at time HI we are describing an exact observable 
on A[- n -i,n+i]- We write (tro-R)^ := Y^=o Rs+4j,s+4k for tracing out the least significant 
qudit, and E 00 = |0)(0| is an elementary 4x4 projector. Note that T n _ m is just a matrix 
of T[ m n ] in the canonical basis \s). It is rather trivial to exactly solve this dynamics for 
a small finite n — m, however this does not yield physically very useful information 
about the KI dynamics. One would wish to study the correct TL by first taking 
n — > oo, and only after that m, t — > oo, however this task seems almost computationally 
intractable. Still, we were able to make some modest numerical experiments exploring 
this question, suggesting that for sufficiently strong integrability breaking (say QC case) 
the asymptotic matrix i?( m '°°)(oo) has a remarkable scale invariance if we coarse-grain 
it by tracing over its 4 x 4 blocks: 

R( m+1 >°°\oo) = tro/^'^oo) = 0R {m,oo) (oo), (42) 

where ( is some scaling factor. This seems to be true for both orders of the limits t — > oo, 
n — ► oo, although better numerical results have been obtained for the 'incorrect' limit, 
namely letting the number of iterations t —>■ oo for a finite register size n, and then 
checking the convergence of results with increasing n. 

Before discussing numerical results we note another useful observation. Let us define 
and briefly study KI chain on a semi-infinite lattice Z + = [0, oo], with the Hamiltonian 
(P]) for L — oo and open boundary condition on the left edge. Now we consider a quasi- 
local algebra A% + . The time automorphism is again strictly local T + : A[q iU ] ~^ «4[o,n+i]> 
and can be written as a semi- infinite product T + = n7ez + Wjj+i- In the definition of 
the truncated time map Tt , = P[ 0jn ]T + : A[o >n ] — > A[o, n ] truncation is needed only on 
the right edge. Note that due to this property, simulation of local observables (localized 
near the edge of the lattice) is twice as efficient than on doubly-infinite lattice, meaning 
that with the same size of computer register one can exactly simulate for twice longer 
times. For example, computation of correlation functions Cg A (t) := (-B|[T + ]M) is exact 

(B\[T + YA) = (B|[T+ n] ]U), for t < 2(n- q), if A, Be A m . (43) 

Numerically inspecting the correlation functions of a simple local spin A = B = ctq 
in fig. [18] we find clear asymptotic exponential decay for QC case, whereas for IN and 

| When labeling tensor product matrix elements we shall always follow a convention that left factors 
are labelled with less significant digits, namely (A <g> B)j + j>d,k+k'd = Aj^By ^> ■ 
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Figure 18. Decay of correlations C + (t) = (cq^ctq) for semi-infinite KI lattice, with 
different cutoff sizes n = 15, 13, 11, and for different cases (QC, NE, IN), all indicated 
in the figure. Note that the two curves for NE case are practically overlapping, and 
that all curves for different r's are exactly overlapping until t — 2r sma n er . 



NE cases we find non-vanishing plateaus in the correlation function, i.e. non-vanishing 
stiffness D + := C + ^ which signals non-ergodicity and existence of local (for integrable 
cases) and pseudo-local (for non-ergodic and non-integrable cases, e.g. NE) conservation 
laws of T + . We note that the asymptotic correlation decay in non-integrable cases, say 
QC and NE, seems quite insensitive to increasing truncation order n - indicating that the 
leading eigenvalues of , remain frozen when increasing n. Note that the asymptotic 
exponents of correlation decay for a semi-infinite chain are not the same as for an infinite 
one, i.e. the point spectra of n , and T n are in general different, however we have some 
indications to believe that their phase diagrams should agree, namely ergodic regimes 
in the KI model on semi-infinite chain model are in one-to-one correspondence with 
ergodic regimes of the model on infinite chain. 

However, the most remarkable feature of dynamics T + is the following. Splitting 
the truncated semi-lattice as [0, n] = [0, m — 1] U [m, n] and following the time evolution 
of an observable A(t) = [Tj[j n i]*A(0) = J2 S a (li s n ) l s °' • • • > s «) * n terms of a 'super-wave- 
function' O/*^ s %, one can again define the reduced super-density matrix as 

R [ r n) u. n(t)= Y af ] ,a?* , , v (44) 

(s m ,...s n ),(s' m ,...,s' n )\ > / j (so,..,s m _i,s m ,...s„) {so,...,s m -i,s' m ,...s' n j v ) 

SQ,---,S m -l 

The dynamical equation for R( m+1 ' n+1 \t + 1) in terms of R^ m ' n \t) is exactly the same as 
for doubly-infinite lattice, namely eqs. fl4QII4ip . Hence also the conjecture (j42p on scaling 
mvanance 

of i?Koo)(oo) should be the same for the two lattice topologies. However, 
numerical computations are much easier and thus the results are more suggestive for 
the semi-infinite case. 
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(a) (b) (c) 




123456 123456 123456 



P P P 

Figure 19. Scaling of the partial norms u(p) = i?Q™ 0+p '™^ (t) in log scale (units 
are arbitrary) versus the partial tracing index p, for cases QC (a), NE (b) and IN 
(c). Diamonds, stars, and squares represent data for the semi-infinite KI lattice with 
truncation size n = 15, (diamonds, stars) and n = 11 (squares), all for too = 4. Finite 
number of time steps t = t* = 17 (diamonds), just before the absorbing boundary 
affects any of the data shown, is compared to steady state observable t — oo (stars). 
In the non-integrable cases (a,b), data are compared also with steady state t = oo 
simulation of the same local initial observable ^4(0) = erg on the two-sided (doubly- 
infinite) KI lattice with truncation size n — 7 (r = 15), and toq = 2 (triangles). 



Let us now discuss some numerical results. We have always started from local initial 
operator A(0) = ctq. We took n as large as allowed by existing computing resources, 
namely n = 15 for the semi-infinite chain and n = 7 for the infinite chain, and compared 
the data for asymptotic matrices i?( m ' n )(oo) (in numerics t has been chosen such that 
the results converged, typically t « 100) with finite time data R^ m,n \t*) where time 
t* was set as large as allowed so that the data were still exact and no truncation was 
needed, typically t* ~ n. In all cases, numerical results were quite insensitive to small 
changes in truncation order n. First we have computed the scaling of the principal 
matrix element, or the partial norms (t) = „ „ ^ \a~p nn J 2 

' ? UU V J ^...Sm-2,Sm-l£a4 1 (. . . ,S m -2 ,S«i- 1 ,0,0. . .) I 

which, assuming (H2l) . should asymptotically scale as oc ( m (see figHHJ). If asymptotic 
dynamics t — > oo is determined by normalizable eigenvectors of T + , which necessarily 
correspond to uni-modular eigenvalues, then we should have ( = 1. In QC case a clear 
scaling was observed for both topologies (Z and Z + ) with the same exponent (, however 
the exponent was slightly different for #( m ' n )(oo) and R^ m ' n \t*). In the cases of non- 
ergodic dynamics (NE and IN) the results for two topologies were quite different. For 
semi-infinite topology we find very clearly that ( = 1 indicating that A(t*) and even 
A(oo) can be written as I 2 convergent sums of local operators. 

As a more quantitative test of conjecture (|42|) we compare the upper-left ('most 
important') 16 x 16 block of the super density matrix scaled to a unit principal element 
Rl k (t) = R { ™ 0+P ' n \t)/R { ™ 0+P ' n) (t). In figj20]we plot the diagonal elements R^/t) for 
different p, while in figj21]we plot 2D charts of the entire scaled density matrices R P k {t). 
Indeed we find for QC case that the matrix R P k (t) is practically insensitive to increasing 
truncation (jp) and to topology of the lattice (semi- infinite versus infinite), both for finite 
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Figure 20. Scaling of the diagonal elements of scaled reduced super-density matrices 
R? j(t) for initial tracing index mo = 6, and p = 1 (diamonds), p = 2 (stars), p = 3 
(squares), p — 4 (triangles), for non-integrable cases of KI chain on semi- infinite lattice 
truncated at n = 15, namely for case QC (a,b), and case NE (d,e). Finite time data 
at t = t* = 17 kicks are shown in (a,d), while steady state observables t — oo are 
analyzed in (b,e), all starting from initial observable CTq. For comparison, asymptotic 
steady state t = oo data for two-sided, doubly infinite KI chain, truncated at r = 13, 
and with initial tracing index mo = 2, are shown in (c,f), namely for QC case (c) and 
NE case (f). Note that data in plot (b) are practically exactly overlapping. 



time t = t* and 'steady-state' t = oo. On the other hand, in non-ergodic cases, the 
scaling ( H2l) is typically broken. However, it seems to be observed in the steady state 
(t = t*) of NE case, which is (in our setting) probably a non-physical but still quite 
robust effect due to truncation (a kind of absorbing boundary condition). 

Summarizing this subsection, we conjectured that in the regime of quantum chaos 
reduced super-density matrices of time-evolved observables typically obey the scaling 
law (H2l) with the exponent ( which only depends on global dynamics on quasi-local 
algebra of observables and not on a particular choice of initial observable. We expect 
that accurate numerical calculations of exponent ( would be possible within a certain 
quantum dynamic renormalization group scheme, however its precise formulation at 
present remains an open problem. 
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Figure 21. Scaling of scaled reduced super density matrices R p - k (t) - the grayness 
level is proportional to \R P k (t)\ - for initial tracing index mo = 4, and p = 1, ... ,6 
(columns plots), for six different cases, QC, NE, and IN of truncated semi-infinite KI 
lattice (n = 15), at finite time t = t* = 17, and asymptotic steady state t = oo (row 
plots), always starting form initial observable (Jg. 16 x 16 'most important' matrix 
elements are plotted with the matrix site j = k = at lower-left corner. 
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7.6. Dynamical entropies 

In section [5J further elaborated in subsection 17.41 we have proposed an entanglement in 
operator space of quasi-local algebra as a possible new measure of quantum algorithmic 
complexity. Here we would like to compare this briefly to some established proposals of 
quantum dynamical entropy, such as for example CNT entropy [19] or AF entropy |20| . 
both being ideally suited for a quantum dynamical system formulated in terms of time 
automorphism T, tracial invariant state u, and quasi-local C* algebra Az- We shall 
here briefly review AF entropy which is conceptually simpler. 

Let us start by taking a set of say k elements A a , a — 1, ... k of A.% which form 
a partition of unity, namely Yl a A* a A a = 1. Following the dynamics up to integer 
time t, A a (t) = T l A a , we construct a set of k l elements depending on a multi-index 
a = (a , • • • ,a t -i), namely 

A® = A ao (0)A ai (l) ■ --A^it - 1) = A ao TA ai T ■ ■ ■ TA at l . (45) 

From the homomorphism property and unitality of dynamics it follows that A$ is also 
a partition of unity, l* 7 ^ = 1- Hence using an invariant state lu one can form 

a positive, Hermitian, trace-one, k l x k l matrix 

P%_ = "([AfrAf). (46) 

pW can clearly be interpreted as a density matrix pertaining to dynamically generated 
partition, and its Von Neumann entropy generated per unit time defines the AF entropy 

S AF = sup lim sup — tr [p {t) log p {t) ] , (47) 

{A a } t^oo t 

where, strictly speaking, supremum over all possible generating partitions has to be 
taken. It is not surprising that practical evaluation of 5*af is impossible except for rather 
trivial cases, such as dynamics generated by shift automorphism Si [20J. However, one 
can easily show that a related linear AF entropy 

Sf = sup lim sup — logtr[pW] 2 , (48) 

{A a } t^oo t 

is tractable much more easily, while the behaviour of S AF and S AF is likely to be similar 
in practice. The key observation is to write the purity P AF (t) = tr[p^] 2 = Yl a a PpjaPoip 
in terms of dynamics over a product algebra Az = Az x Az, with time automorphism 
T(AxB) = T(A) x T(B) and an invariant state Co(AxB) = u{A)uj(B). To this product 
structure we add a linear map K : Az — > Az depending on a generating partition {A a } 

k 

K(B) = i£fi x A «)B{K x A* p ). (49) 

o,/3=l 

If 1 is a unit element in A% then 1 x 1 is a unit element in Az, and purity can be written 
using a transfer-matrix-like approach as 

P AF (t) = Cj |[TK]*(1 x 1)| . (50) 
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The idea can be worked out in detail for a complete generating partition of a local sub- 
algebra A[-. q +i >q ] of size k = 4 2q . There it turns out that, due to locality of dynamics 
(04"]) . the resulting purity is independent of q > 1, i.e. the supremum f T4Hl) is already 
achieved by a rather modest partition of 4 2 = 16 elements, and the map TK can be 
factored into a direct product of two maps acting separately on two independent copies of 
a semi-infinite lattice. Leaving out some technical details of derivation, the final result 
reads as follows. Let us consider a time dependent 4* x 4* matrix Q(t), representing a 
state on Az+ = Az+ x A%+, with initial value Q(0) = Qo,o(0) = 1, and dynamics given 
by the following completely positive matrix map 

Q(t + 1) = tr [0)1] {T t+1 [E Q0 <g> 1 4 ® Q(t) <g> E 00 ] I^i} , (51) 

where (tr^iii?)^ = ^ s =o -^s+i6j,s+i6£: traces out two least important qudits, and unitary 
time evolution matrix T t is given in ( T4TT) . Then the purity, and linear dynamical entropy 
(LDE), are simply given as P AF (t) = [Qo,o(t)} 2 , S 2 (t) = — 2 log<5o,o(^)> respectively. The 
asymptotic increase of ^(t) per unit time yields the linear AF entropy ( T4B1) . Again, 
for practical calculations it is convenient to consider truncation of matrices Q(t) after 
each iteration (15TT) . say at dimension 4 r . In fig. [22] we plot LDE for different cases of 
KI dynamics, and we observe that LDE is always clearly growing linearly oc t, so the 
AF entropy is always strictly positive, even in non-ergodic (NE) and integrable (IN) 
cases, and that the results are robust and stable against changing the truncation order 
r. Perhaps this finding appears surprising, but one has to bear in mind that AF and 
CNT entropies can be positive even for simpler dynamics, such as quasi-free flows on 
C* algebras. 

For ergodic dynamical systems on C* algebras one can use Shannon-Mcmillan- 
Breiman (SMB) theorem (for classical SMB theorem see e.g. [78J, and [79J for a possible 
quantum extension), which states that for typical sequences a, multi-time correlation 
function (MTCF) should decay exponentially 

C(t) = v(A®) ~ exp(-ht) (52) 

where the exponent h, which should be essentially independent of a, is equal to 
a dynamical entropy of the map T with respect to an invariant state uj. For an 
interesting application of SMB theorem in the context of quantum dynamical chaos 
see Ref. [80]. In our numerical experiments we considered truncated dynamics T [_ n ,n]> 
writing truncation order as r = 2n + 1, and computed two kinds of MTCF: (i) For 
a uniform sequence a = (1,1,...), where A 1 = Uq (or A\ = o-qO~* in which case 
the truncated lattice was placed as [—n,n + 1], hence r = In + 2), we estimated 
MTCF fl52]) as (^4i|Tr_ n) „]AiT[_ n)n i • • - T\ n ^Ax). (ii) For a random sequence a, we 
took a complete generating partition . . . s q )} on A[- qiq ] with k = 4 l elements, for 

I = 2q + 1 = 1 or I = 3, and computed an average square modulus of MTCF (\C(t)\ 2 ) = 
(|(A a jT[_ n>n ]A ai T[_ rV j] • • • T[_ n<n ]A at _J| 2 ) over 20 randomly sampled sequences a. In 
both cases (i,ii) we have found a rather good agreement between log|C(t)| 2 and LDE 
S*2(t) for the case of quantum chaotic dynamics (fig. 1221) . For average random MTCF we 
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t 



Figure 22. Linear dynamical entropies (continuous curves) for cutoff orders at 
r = 6, 7, 8, for three different cases of dynamics (QC (a), NE (b) and IN (c)), compared 
to - In |C(t)| 2 where C(t) are MTCF with obervables X = of and XX = u\u\ at 
different cutoff orders r (symbols), and to average — ln(|C(£)| 2 ) over random MTCF 
sampled over 20 sequences of random Pauli observables of lengths / = 1 and I = 3 
(curve-symbols) . 
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found rather good agreement even for non-ergodic cases (NE and IN), however MTCF 
for the uniform sequence exhibited big oscillating fluctuations there, indicating simply 
that SMB theorem does hold for non-ergodic dynamics of KI model. 

Summarizing, it seems that AF entropy, even though it very cleanly generalizes 
the concept of Kolmogorov- Sinai classical dynamical entropy to quantum dynamical 
systems, cannot be used as an indicator of quantum chaos or an indicator of 
computational complexity of quantum dynamical systems. Nevertheless, it has to be 
stressed that quantum dynamical entropies can be related to the notion of quantum 
algorithmic complexity, mentioned in section [51 by a quantum version of Brudno theorem 
established in Ref. [81 J . 

8. Conclusions and open questions 

In the present article we have reviewed several possible approaches to describe dynamic 
instabilities, relaxation phenomena, and computation complexity in the simple model 
of one-dimensional non-integrable locally interacting quantum many body dynamics. 
We have argued that the essential non-equilibrium statistical properties of quantum 
dynamical systems - in the absence of idealized external baths - may be crucially related 
to the integrability of the system, or in the complementary case, to the existence of the 
regime of quantum chaos. 

The general flavor which remains after such studies is that the non-integrable 
quantum many body problem at high temperature will preclude any exact and complete 
solution by its very nature. Still it is hoped that a more complete ergodic theory of such 
systems could be developed, allowing for example for exact calculation of relaxation 
rates, scaling exponents of resonance eigenvectors, etc. 

It has been shown in many cases that integrable quantum systems have anomalous 
non-equilibrium statistical mechanics at high temperature, for example they exhibit 
ballistic transport. This can easily be understood as being a consequence of existence 
of (an infinite sequence of) exact local conservation laws which prevents quantum 
ergodicity, similarly as existence of canonical action variables in classical integrable 
systems prevents classical ergodicity. An interesting open question is the following: 
How strong integrability breaking perturbation is needed, in generic cases, to break all 
the exact conservation laws and yield normal (diffusive) transport? In other words: A 
quantum KAM theory is needed! Numerical experiments shown in this paper suggest an 
interesting possibility, namely that in some cases a finite, critical perturbation strength 
is required. 
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